Claude Code Review

Head SHA: 9879df3

Files changed:

• 21
• src/simulation/m_ibm.fpp
• src/simulation/m_start_up.fpp
• src/simulation/m_collisions.fpp
• src/simulation/m_global_parameters.fpp
• src/common/m_constants.fpp
• src/common/m_derived_types.fpp
• src/common/m_mpi_common.fpp
• src/post_process/m_data_output.fpp
• src/simulation/m_ib_patches.fpp
• src/simulation/m_time_steppers.fpp

Findings:

Critical: Stack overflow in s_reduce_ib_patch_array (src/simulation/m_start_up.fpp)

subroutine s_reduce_ib_patch_array()
    type(ib_patch_parameters), dimension(num_ib_patches_max) :: patch_ib_gbl

patch_ib_gbl is a local (stack-allocated) variable dimensioned by num_ib_patches_max, which was just raised from 50,000 to 2,050,000 in src/common/m_constants.fpp. ib_patch_parameters contains a character(LEN=pathlen_max) field (pathlen_max=400) plus ~500 bytes of numeric fields — roughly 900 bytes per element. At 2,050,000 elements this is ~1.8 GB on the call stack, which will segfault on every platform (typical stack limits are 8–64 MB). This variable must be heap-allocated (allocate/deallocate).

High: Debug print left in production code (src/simulation/m_ibm.fpp)

print *, proc_rank, " New Owner ", patch_ib(k)%gbl_patch_id  ! TODO :: REMOVE THIS DEBUG PRINT

This fires on every ownership transfer for every locally-tracked IB, generating unbounded output at scale. The ! TODO :: REMOVE THIS DEBUG PRINT marker confirms it is unintentional.

Medium: Commented-out code in m_collisions.fpp leaves pid2 lookup absent (src/simulation/m_collisions.fpp)

! call s_get_neighborhood_idx(pid1, pid1) ! global patch ID -> local index call s_get_neighborhood_idx(pid2, pid2)
if (pid1 <= 0 .or. pid2 <= 0) cycle

The comment text contains call s_get_neighborhood_idx(pid2, pid2) — what appears to be a second intended lookup that was accidentally folded into the comment instead of being left as executable code. Neither lookup actually executes: pid1 and pid2 are the raw decoded global IDs from s_decode_patch_periodicity, not local indices. The guard if (pid1 <= 0 .or. pid2 <= 0) cycle would never trigger for valid global IDs (which are ≥ 1), meaning the subsequent patch_ib(pid1) and patch_ib(pid2) accesses use global IDs as local array indices, which is out-of-bounds when num_ibs < num_gbl_ibs. If this is intentional (global IDs equal local indices in the no-MPI / single-rank path), that should be documented; otherwise both lookups need to be uncommented.

Low: GPU_UPDATE(device=[patch_ib(ib_idx)%moment]) removed without replacement (src/simulation/m_ibm.fpp)

-            patch_ib(ib_marker)%moment = moment*patch_ib(ib_marker)%mass/(count*cell_volume)
-            $:GPU_UPDATE(device='[patch_ib(ib_marker)%moment]')
+            patch_ib(ib_idx)%moment = moment*patch_ib(ib_idx)%mass/(count*cell_volume)

The per-field device update was removed. If patch_ib(i)%moment is consumed on the GPU before the next full GPU_UPDATE(device=[patch_ib(1:num_ibs)]) (e.g., in a subsequent call to s_ibm_correct_state within the same time step), the device will use a stale value. Verify that a covering full-struct update always precedes any GPU use of %moment after s_compute_moment_of_inertia is called from s_propagate_immersed_boundaries.

You had an issue that I think I fixed.

Root cause: num_ib_patches_max = 2,050,000 caused a ~2.25 GB unconditional allocation + full init loop on every startup (even for non-IBM cases), crashing CI debug builds.

Fix — 3 files, 2 commits:

1. m_constants.fpp: Added num_ib_patches_max_namelist = 50000 — the small initial budget (same as before particle beds were added).
2. m_global_parameters.fpp + m_start_up.fpp: Initial patch_ib allocation uses num_ib_patches_max_namelist on all ranks (~55 MB instead of ~2.25 GB).
3. m_particle_bed.fpp: Grows patch_ib to num_ib_patches_max via MOVE_ALLOC only when particle beds are actually being generated (rank 0 only, before the MPI broadcast).
4. m_mpi_proxy.fpp: Non-rank-0 processes grow patch_ib to match if num_ibs > size(patch_ib) — handles the MPI > 1 case where rank 0 grew for particle beds but other ranks still
   have the small allocation.
5. m_start_up.fpp s_reduce_ib_patch_array: Changed the 2.25 GB stack-allocated patch_ib_gbl to a heap-allocated array sized to exactly num_ibs — eliminates the SIGSEGV for IBM
   cases.

I generated and pushed the missing golden file for BA186DDF (2D -> Example -> mibm_particle_bed) directly to your branch — that should fix the CI failure.

While reviewing the diff I found several other issues worth addressing:

Critical (wrong results)

f_neighborhood_ranks_own_location guard is num_procs > 2 instead of num_procs > 1 (m_collisions.fpp)
On a 2-rank run the condition is false, so both ranks think they own every collision and forces get double-applied. Should match f_local_rank_owns_location's guard.

recv_forces_snap not zeroed before y/z accumulation passes (m_ibm.fpp, s_communicate_ib_forces)
recv_forces_snap = 0._wp lives only inside the x-dimension block. For y and z it accumulates stale values from the previous dimension, making the double-count subtraction remove too much force.

s_handoff_ib_ownership can silently drop newly received patches (m_ibm.fpp)
After the compact-and-delete loop, ib_gbl_idx_lookup still reflects the old layout (it is rebuilt after unpacking). s_get_neighborhood_idx is therefore testing the pre-handoff table — a just-evicted patch may still appear as resident, causing the newly received authoritative state to be discarded.

High (MPI / GPU correctness)

recv_forces_snap, recv_torques_snap, recv_ids, recv_ft allocated with plain allocate, not @:ALLOCATE (m_ibm.fpp)
These arrays are used inside GPU_PARALLEL_LOOP via copyin clauses, but plain allocate skips GPU_ENTER_DATA. On OpenACC builds the GPU sees uninitialized device memory.

patch_ib re-allocated with plain allocate/deallocate while GPU_DECLARE is active (m_global_parameters.fpp + m_start_up.fpp)
patch_ib was changed from a static array to allocatable, but the device pointer registered by GPU_DECLARE at startup goes stale every time s_reduce_ib_patch_array or the MPI proxy resizes it. Each reallocation should use @:ALLOCATE/@:DEALLOCATE.

neighbor_domain_x/y/z computed on host but never GPU_UPDATE'd to device (m_start_up.fpp)
These are GPU_DECLARE'd and used inside f_neighborhood_ranks_own_location (a GPU_ROUTINE). After get_neighbor_bounds fills them on the CPU, a GPU_UPDATE(device=[neighbor_domain_x, neighbor_domain_y, neighbor_domain_z]) is needed.

s_compute_ib_neighbor_ranks: nreqs not reset between outer loop iterations (m_start_up.fpp)
After MPI_WAITALL, nreqs retains its end value going into the next k iteration. With ib_neighborhood_radius > 2 this overflows the dimension(52) requests array.

Post-process SILO writes called from all ranks (post_process/m_data_output.fpp)
DBPUTMVAR to dbroot is now called from all ranks, not just rank 0. In file_per_process=F mode this produces duplicate or corrupted multi-mesh entries.

Medium (robustness / performance)

num_local_ibs_max = 2000 hard cap with no overflow check — local_ib_patch_ids(num_local_ibs) = num_ibs writes out of bounds silently if more than 2000 particles land on one rank. Needs an @:PROHIBIT.

ib_force_send_buf/recv_buf heap-allocated on every RK stage call (m_ibm.fpp) — these should be module-level buffers sized once in s_ibm_setup rather than allocated and freed on every timestep.

moving_immersed_boundary_flag not recomputed after s_handoff_ib_ownership — set once at startup, then stale after each ownership handoff cycle when moving IBs migrate between ranks.

Low (code quality)

• s_get_neighborhood_idx has an unused local variable i (compiler warning)
• f_xorshift can return exactly 1.0 (violates the [0, 1) contract in the docstring — particle placed at xmax)
• particle_bed_parameters type defined in simulation-only m_global_parameters.fpp; convention is to put derived types in m_derived_types.fpp

The two highest-priority fixes are the num_procs > 2 off-by-one (one-character fix, guaranteed wrong on 2-rank jobs) and recv_forces_snap not being zeroed per-dimension (wrong forces in any 2D/3D multi-rank run).

Addressing AI comments:

All of the critical issues are incorrect for one reason or another. I have gone through and verified that the relevant chunks of code are in place to address the concerns it has.

As for the neighbor boundaries, they are never called in a GPU region, so this comment is mute. I have removed the GPU routine call to help with compilation in response.

The file per process F claim is just false, and I tested that the code runs find without it. But it is true that we do not need other ranks to be doing any of this calculation regardless. I extended the rank 0 check in response.

nreqs is in fact reset between loop iterations.

Snapshot arrays are in fact reset between loop iterations

It is confused about the need to GPU allocate arrays for patch_ib. Because the array is used with a declare, but never actually allocated on the GPU, we just have the GPU pointer. When we perform the first copy to the GPU to allocate the array size, we have already done the reshape. This means that the GPU array is properly sized. This is important because the size of the patch_ib array is significant in the new code and could causes issues with limiting problem size because the GPU must be able to hold the global patch_ib array, even though it would then immediately be deallocated and reallocated to a smaller array. If the number of ib patches grows to the extent that it costs many for GB of memory, this could cause memory limitations. The current implementation should be the correct way to go.

the moving IBM flag should not be recomputed. You want it to be the same on all ranks. There is argument that it should be moved to the case level for case optimizaiton later, but that is a future PR.

Force arrays getting moved is going to be handled in a future PR that already has a github issue I created yesterday.

Prohibit for number of local IBs resolved.

In summary, two of these I made minor adjustments for. One I reject for reasons that I believe are valid and better for long-term code support. All others are incorrect.

General Note: Mac errors are on sections of the code that were not touched by this PR (bubbles), and it is a random dyl linker error. I get the same error on reldebug on a chemsitry case. I get the same error locally, but on a different set of tests. This seems like a MacOS instability issue, and not related to the particular PR. Going to attempt to rerun. If ti fails again, but on different tests, then it looks like the MacOS test suite checks may be broken.

Edit: Rerunning the job changed the jobs which failed. This appears to be a spurious failure for MacOS specifically. It appears to be specific to the branch, as other branch are not experiencing it. Will investigate.